TURBULENCE  MODELING  VALIDATION  FOR  AFTERBODY  FLOWS 
A.  Hadjadj1  and  A.N.  Kudryavtsev2 

1  LMFN-CORIA,  UMR  6614,  INSA  de  Rouen,  B.P.  8, 

76801  Saint-Etienne  du  Rouvray,  France 
2  Institute  of  Theoretical  and  Applied  Mechanics, 

630090  Novosibirsk,  Russia 

Introduction 

The  aerothermodynamics  phenomena  over  space  launch  vehicles  or  missiles  are  a 
challenging  problem  on  space  and  aeronautical  applications.  These  physical  phenomena  can 
strongly  affect  the  engine's  aerodynamic  performances.  The  physical  problem  met  on  this 
geometry  is  essentially  the  result  of  the  interaction  of  two  merging  flows;  one  issued  from  a 
propulsive  jet  at  high  speed  and  high  temperature  (flow  2),  and  the  other  (flow  1 )  caused  by  the 
ambient  low  speed  stream  (see  Fig.  1). 

Several  complicating  issues  are  presented  in  this  flow  configuration  : 

•  The  afterbody  flow- field  is  dominated  by  the  behavior  of 
the  turbulent  boundary  layer.  Thus  the  quality  of 
computational  results  depends  strongly  on  the  accuracy  of 
the  turbulence  model  used  for  CFD. 

•  The  afterbody  configurations  are  often  quite  complex, 
since  they  include  single-  or  double  flux  nozzles  with  a 
variable  geometry,  twin-engine  with  tails  and  non 
axisymmetric  nozzles  with  base  flows.  All  of  these 
elements  are  important,  and  should  be  modeled. 

•  On  the  other  hand,  a  rather  fair  prediction  of  the  flow- 
field  can  be  accompanied  by  large  errors  in  the  calcula¬ 
tion  of  wall  properties,  affecting  mainly  transfer  coeffi¬ 
cients,  skin  friction  and  heat  flux. 

The  complexity  and  the  usual  misunderstanding  of  these  physical  phenomena  coupled 
with  the  turbulence  modeling  problem  of  complex  compressible  flows  motivates  the  present 
study.  The  outline  of  this  paper  is  as  follow:  the  turbulence  models  used  in  this  study  are  given 
in  section  2,  the  main  lines  of  the  numerical  method  are  presented  in  section  3,  the  nozzle  start¬ 
up  process  is  examined  in  section  4,  finally  sections  5  and  6  are  devoted  to  the  comparisons 
between  numerical  and  experimental  measurements.  It  is  found  that  an  interesting  estimation  of 
the  pressure  coefficient  on  the  afterbody  is  obtained  via  the  multi-scale  model 

Turbulence  modeling 

The  oldest  and  simplest  proposal  for  modeling  the  turbulent  Reynolds  stress  is  Boussinesq 
isotropic  eddy-viscosity  concept  that  assumes  an  analogy  between  the  viscous  stresses  in 
laminar  flows  and  the  turbulent  stresses  in  turbulent  flows.  This  concept  relates  linearly  the 
Reynolds  stresses  to  the  mean  flow  quantities  as  : 

p U-Uj  =  2/3  p  k  Sff  -  H '  (, Stj  -  1/3  Sa  5 ,j)  (1) 

where  Sy  represents  the  main  strain  rate  and  p  t  is  the  eddy  viscosity  defined  by  : 

p  t  =  p  k2  /  c .  (2) 


(i) 


Figure  1.  Afterbody  nozzle 
configuration. 
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Equation  (2)  constitutes  a  common  basis  for  most  turbulence  models  extensively  used  today  [6]. 
The  baseline  model  of  this  study  is  the  standard  k-z  model  proposed  by  Jones  and  Launder  [3]. 
This  model  based  on  Boussinesq  eddy  viscosity  concept  (relation  2)  uses  two  supplementary 
equations  to  describe  the  whole  turbulent  flow  field. 

In  this  paper,  we  propose  the  introduction  of  several  scales  characterizing  different 
process  of  turbulence  interaction  by  using  multi-scale  model  [2],  [4],  [7].  The  multi-scale 
turbulence  modeling  consists  in  partitioning  the  energy  spectrum  into  several  regions,  each 
characterized  by  a  different  time  scale.  Thus  different  mechanisms  such  as  the  return  to 
isotropy,  dissipation,  diffusion,  and  viscosity  can  be  modeled  using  the  various  characteristic 
scales.  In  the  simplified  split-spectrum  used  by  Kim  [4],  turbulent  transport  is  described  by 
using  two  time  scales:  the  first  one  corresponds  to  the  large  eddies  and  describes  generation  of 
turbulent  kinetic  energy,  and  the  second  one  corresponds  to  the  smaller-scale  eddies  and 
describes  dissipation  rate.  The  turbulent  energy  can  be  partitioned  into  three  regions:  production 
region,  characterized  by  the  turbulent  kinetic  energy  kp  and  the  energy  transfer  rate  z  p,  transfer 
region,  characterized  by  the  turbulent  kinetic  energy  kt  and  the  dissipation  rate  z  t  and 
dissipation  region,  where  the  turbulent  kinetic  energy  is  dissipated  into  heat. 

This  model  employs  four  supplementary  transport  equations  to  described  large  and  small 
eddies  motion  (for  more  details  see  Kim  [4]).  One  of  the  most  interest  of  the  multi-scale  model 
is  the  expression  of  the  eddy  viscosity,  given  by: 

P  t  =  Q,  P  k2  /  e  p .  (3) 

where  k  =  kp+  kt.  By  rearranging  the  equation  (3),  we  obtain  : 

p  /  =  Cp  (e/e p)  pk2  /  e  ,=  F(e,/ep)  p  k2  /  e  .  (4) 

Note  that,  in  the  proposed  multi-scale  model,  becomes  function  of  zt/zp.  Therefore,  for 
turbulence  equilibrium  conditions,  this  function  F  is  equal  to  a  constant  value  in  the  classical  k- 
e  model.  In  this  case,  the  multi-scale  turbulence  model  becomes  a  single  scale  turbulence 
model.  In  the  case  of  non  equilibrium  turbulent  flow,  for  example  when  the  energy  transfer  rate 
from  large  eddies  is  more  important  than  the  dissipation  rate  at  small  eddies,  the  eddy  viscosity 
p  t  decreases,  which  implies  a  decreasing  of  the  turbulent  kinetic  production.  Inversely,  this 
term  permits  to  increase  the  turbulent  kinetic  production.  The  multi-scale  spectral 
decomposition  yields  to  a  variable  Cp  which  accounts  for  turbulence  non  equilibrium  effects. 

Computational  results 

Afterbody  mesh  grid.  The  computational  domain  is  splitted  into  two  connected  sub- 
domains  (see  Fig.  2)  :  the  first  one  includes  the  nozzle  with  a  propulsive  jet,  and  the  second  one 
concerns  the  ambient  low  speed  flow.  The  downstream  region  of  the  nozzle  has  been  described 


Figure  2.  Computational  domain  (Grid  detail  near  trailing  edge). 
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in  detail  with  a  sufficient  number  of  points  in  order  to  solve  correctly  the  throat  region  where 
sonic  speed  is  attained.  Additional  mesh  grid  refinement  was  provided  near  the  trailing  edge  of 
the  nozzle  so  as  to  predict  as  well  as  possible  the  entrainment  effect  issued  from  boundary 
layers  evolution,  into  the  free  shear  layer  between  the  plume  and  the  free  stream  flow.  The 
computational  domain  is  discretized  into  a  330x150  mesh  with  a  concentration  of  grid  points 
near  afterbody  walls.  The  mesh  size  of  the  first  grid  point  from  the  wall  is  located  at  y+ &  1. 
Extensive  grid  refinement  tests  were  achieved  to  insure  that  the  results  are  mesh  independent. 

Numerical  aspects .  The  Reynolds-averaged  Navier  -  Stokes  equations  associated  with  a 
turbulence  model  are  solved  on  a  computational  domain  of  variables  £,  and  r|  (transformed  coor¬ 
dinates  of  the  physical  domain),  by  the  use  of  finite  volumes  discretization  technique  on  struc¬ 
tured  mesh.  The  new  system  of  equations  is  solved  by  using  an  explicit-implicit  scheme  which 
is  second-order  accurate  in  space  and  time  (see  Vandromme  et  al.  [8],  [9]).  The  basic  discretiza¬ 
tion  for  the  convective  fluxes  is  modified  so  as  to  take  into  account  the  physical  properties  of 
information  propagation.  The  viscous  terms  are  centered  and  the  axisymmetric  source  terms  are 
integrated  at  the  center  of  each  control  volume  in  both  the  £,  and  r|  directional  sweeps. 

To  reach  a  steady-state  solution  with  a  minimum  number  of  iterations,  the  explicit 
discretization  is  complemented  with  an  implicit  numerical  approximation  which  is  free  from 
stability  conditions.  Thus,  the  block  pentadiagonal  system  is  solved  by  generalized  Thomas 
algorithm  with  LU  decomposition  in  the  r|  direction,  and  by  a  line  Gauss-Seidel  relaxation 
technique  in  the  £,  direction.  As  the  system  is  really  diagonally  dominant  and  the  method  is 
iterative,  converged  steady  solution  can  be  obtained  in  very  limited  number  of  time  steps,  each 
time  step  including  a  double  sweep  (backward-forward)  in  the  flow  direction.  With  this 
technique,  unbounded  time  step  values  can  be  used.  Numerical  simulations  have  been  made 
with  C.F.L.  numbers  greater  than  104. 

The  numerical  boundary  conditions  used  in  this  work  are  far-field  non-reflecting 
boundary  conditions  for  ambient  low  speed  stream,  adiabatic  no-slip  wall  and  symmetry 
boundary  condition.  At  the  inlet,  the  profiles  of  all  quantities  are  imposed;  including  total 
pressure  and  temperature  Ph  7),  turbulent  kinetic  energy  k  and  dissipation  rate  e.  Theses  profiles 
are  based  on  the  computational  results  of  the  fully  developed  boundary  layer  flow. 

Axisymmetric  afterbody  flow  (hot  jet).  In  this  study,  we  are  interested  in  the  S3Ch  single¬ 
flux  axisymmetric  nozzle  [1],  [10],  for  which  detailed  measurements  are  available.  This  test 
case  permits  an  accurate  validation  of  turbulence  models  using  CFD.  The  nominal  flow 
conditions  are  : 

•  Outer  flow 

-  upstream  Mach  number  Maj  =  0.8 

-  stagnation  pressure  :  Pt  =  105  Pa 

-  stagnation  temperature  :  Tt  =  300  K 

-  incident  boundary  layer  thickness  (5^8  mm) 

In  order  to  analyze  the  basic  characters  of  the  physical  phenomenon,  both  steady  and  unsteady 
calculations  for  inviscid  and  viscous  flow  have  been  performed. 

Nozzle  start-up  process 

Firstly,  it  is  interesting  to  analyze,  via  numerical  simulation,  the  nozzle  start-up  process. 
This  simulation  starts  when  the  diaphragm,  located  at  the  inlet  nozzle,  is  removed  till  the 
establishment  of  the  propulsive  jet  far  field  downstream  of  the  nozzle.  The  calculation  is 
performed  by  solving  Euler  equations,  since  the  inviscid  effect  (waves  propagation  and  vortex 
formation)  plays  a  major  role  during  this  phase. 


•  Jet  flow 

-  stagnation  pressure  :Pt  =  3.15  105  Pa 

-  stagnation  temperature  :  Tt  =  900  K 
-jet  expansion  ratio  :  Pt!Pa  =  4.80 
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Figure  3.  Pressure  field  (min  =  0.65,  max  =  3.70)  bar. 


Figure  4.  Pressure  field  at  t  =  0.20  ms 
(min  =  0.10,  max  =  3.59)  bar 


RECOMPRESSION  SHOCK 

Figure  5.  Pressure  field  at  t  =  0.35  ms 
(min  =  0.09,  max  =  4.21)  bar 
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As  mentioned  above,  the  computation  starts  impulsively,  and  the  pressure  ratio  is  assigned 
to  be  4.80.  Several  time  sequences  of  color  pictures  (pressure  field)  are  presented  in  Figs.  3-5  to 
illustrate  the  shock  waves  propagation.  Figure  3. a  depicts  the  starting  shock  pulsed  within  the 
nozzle.  Under  a  high  pressure  gradient  effect,  and  due  to  nozzle  geometry  (convergent  part)  the 
incident  shock  interacts  with  the  wall  and  creates  a  series  of  complicated  waves  ;  A  reflected 
shock  propagates  upstream  of  the  nozzle,  a  quasi-normal  shock  near  the  wall  and  a  slip  line 
emanates  from  the  triple  point  (see  Fig.  3  .b). 

The  new  shock  configuration  interacts  with  the  axis  and  generates  a  complex  network  of 
waves  (see  Figs.  3.c  and  4).  At  the  nozzle  exit,  the  flow  is  highly  under-expanded,  and  conse¬ 
quently,  a  low  density  vortex  core  is  generated  (see  Fig.  4).  The  reversed  flow  is  afterwards 
convected  downstream  during  the  starting  process.  Upstream  of  the  incident  shock,  the  flow  is 
highly  decelerated  and  the  static  pressure  is  nearly  equal  to  the  stagnation  pressure.  In  the  super¬ 
sonic  jet,  a  double  compression  process  is  achieved  through  a  straight  longitudinal  Mach  disc 
and  a  curved  transversal  shock.  At  this  time,  we  observe  once  again  a  brief  apparition  of  the  tri¬ 
ple  point  configuration  (see  Fig.  5).  In  fact,  this  flow  configuration  is  unstable.  Thus,  at  t  =  0.3 
ms,  the  Mach  disc  disappears  and  the  oblique  shock  deviates  towards  the  axis.  The  inclination 
of  the  stream  causes  the  shock  wave  to  be  deflected  and  meet  the  axis  at  a  single  point. 

Finally,  two  different  streams  are  formed;  the  first  one  represents  the  propulsive  jet  with 
downstream  alternatives  expansion  and  compression  regions  and  the  second  one  constitutes  the 
uniform  subsonic  flow. 

In  the  following  sections,  we  focus  our  study  on  the  turbulence  modeling,  our  main 
objective  is  to  seek  out  the  possibilities  of  using  the  multi-scale  turbulence  model  to  reproduce 
the  major  features  of  nozzle  afterbody  flows. 

Flow  feature  description 

For  a  laminar  flow,  the  velocity  field  shows  a  separated  boundary  layer  along  the  nozzle 
shoulder.  This  result  is  illustrated  by  the  velocity  vector  field  near  the  nozzle  trailing  edge  (see 
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Figure  6.  Velocity  vector  on  the  afterbody/ 


Fig.  6.  a).  The  boundary  layer  separation  phenomenon  has  been  already  observed  in  many 
calculations  mainly  those  done  by  SNECMA  and  Me  Donnell  Douglas  [10],  [1].  However,  this 
result  do  not  reflect  reality.  Indeed,  the  experimental  boundary  layer  velocity  profiles  indicate 
that  there  is  no  separation  along  the  nozzle  shoulder  [1]. 

The  obtained  results,  even  if  they  are  in  contradiction  with  the  experimental  data,  can  be 
easily  explained  by  the  boundary  layer  effects.  It  is  clear  that  the  laminar  boundary  layer  profile 
developed  along  the  nozzle  shape  is  less  energetic  and  consequently  less  resistant  to  the  adverse 
pressure  gradient,  which  causes  its  separation  from  the  wall.  For  turbulent  flow,  the 
computational  result  shows  that  there  is  no  boundary  layer  separation  on  the  afterbody  (see 
Fig.  6.  b).  This  result  is  in  good  agreement  with  experimental  data. 

Aerodynamics  of  afterbody 

One  important  characteristic  of  predicting  afterbody  flows,  lies  in  the  fact  that  the  surface 
pressure  distribution  may  significantly  affect  the  engine's  aerodynamics  performances.  So,  it  is 
necessary  to  accurately  predict  the  pressure  on  the  afterbody.  Therefore,  the  pressure  coefficient 
Cp  is  given  by  : 

Cp  =  0.5(p-Px)/  pxU\  (5) 

where  the  subscript  "oo"  refers  to  the  upstream  flow  conditions.  The  experimental  afterbody  model 
is  equipped  with  a  total  of  48  pressure  orifices  located  along  the  nozzle  contour.  The  nozzle  and 
the  afterbody  are  defined  in  Fig.  2.  The  geometry  contour  of  the  afterbody  is  made  of  rectilinear 

segments,  the  shape  consisting  of  two  successive  coni¬ 
cal  surfaces  with  an  angle  7°  and  13°,  respectively. 
The  total  length  of  the  afterbody  is  21 8  mm. 

Approaching  the  afterbody  and  nozzle,  the 
boundary  layer  is  growing  slowly  in  a  near  constant 
pressure.  The  shoulder  at  the  beginning  of  the  nozzle 
creates  a  locally  minimum  pressure.  Past  this  mini¬ 
mum  point,  the  gradient  becomes  adverse  and  the 
boundary  layer  begins  thickening  rapidly.  This  adverse 
pressure  gradient  still  acts  until  reaching  the  nozzle 
exit.  The  interaction  in  the  final  adverse  pressure  gra¬ 
dient  region  is  recognized  as  the  key  problem  in  pre¬ 
dicting 

Figure  7  shows  the  calculated  and  the  measured 
wall  pressure  coefficient  Cp  on  the  afterbody.  From  a 


Figure  7.  Static  pressure  coefficient  on  the 
afterbody 
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general  point  of  view,  the  two  simulations  presented  herein  (k-e  and  multi-scale  models),  give 
more  or  less  the  same  trend:  a  rather  good  prediction  of  the  pressure  distribution  in  the  upstream 
part  of  the  afterbody.  However,  at  the  end  of  the  shoulder,  the  k-s  model  shows  an  important  dis¬ 
crepancy  with  the  experimental  data.  The  Cp  is  overestimated  of  almost  30  %,  which  is  unaccept¬ 
able  from  aerodynamics  point  of  view.  This  result  shows  one  of  the  weakness  of  the  k-s  model,  to 
reproduce  the  entrainment  effect  of  boundary  layer  flow  at  the  end  of  the  body.  It  should  be 
pointed  out  that  without  an  accurate  prediction  of  the  pressure  distribution  in  this  region,  any  drag 
estimation  is  meaningless.  Therefore,  the  obtained  results  with  the  multi-scale  model  are  signifi¬ 
cantly  in  better  agreement  with  the  experimental  data.  This  model  reproduces  correctly  the  flow 
behavior  in  this  region;  the  physical  processes  of  boundary  layer  expansion  and  compression  is 
quite  well  described,  and  the  difference  observed  at  the  nozzle  trailing  edge  is  reduced  to  10  %. 

At  this  stage,  it  would  be  important  to  know  that  the  upstream  boundary  condition  have  a 
great  influence  on  the  pressure  distribution  on  the  afterbody.  The  upstream  boundary  of  the 
computational  domain  is  located  in  the  cylindrical  section  of  the  afterbody  model.  The  only  ex¬ 
perimental  data  available  in  this  region  is  a  pitot  pressure  profile.  Without  any  further  informa¬ 
tion  about  the  temperature  and  turbulent  quantities  profiles  in  the  incoming  boundary  layer, 
some  assumptions  should  be  made  to  evaluate  the  missing  quantities.  The  computations  done  by 
SNECMA  [5]  show  that  the  variation  of  the  incident  boundary  layer  modifies  significantly  the 
pressure  value  at  the  nozzle  trailing  edge.  For  example,  a  boundary  layer  with  a  thickness  5  = 
1.5  cm  gives  a  result  closer  to  experimental  point  than  one  with  a  thickness  5  =  0.8  cm.  In  this 
study,  our  strategy  of  investigation  is  to  keep  constant  the  incident  boundary  layer  thickness  6 
and  to  improve  the  accuracy  of  turbulence  modeling. 

Conclusion 

In  this  study,  a  multiple-time-scale  model  was  investigated.  The  basic  form  of  model 
equations  were  introduced  and  applied  to  an  afterbody  nozzle  flow  with  a  propulsive  jet.  It  was 
shown  that  the  present  model  can  satisfactorily  reproduce  the  physical  phenomena  encountered 
in  this  geometry.  The  upstream  conditions  even  if  difficult  to  determine,  have  to  be  considered 
with  accuracy.  This  particular  point  makes  the  computation  even  more  difficult  because  the 
steady  state  depends  strongly  on  the  slightest  variation  of  the  subsonic  boundary  condition. 
Convergent-divergent  axisymmetric  exhaust  nozzles  may  present  considerable  loss  of  thrust  if 
they  operate  under  highly  over-expansion  regime.  Thus  different  nozzle  shapes  like  double-flux 
geometry  will  be  a  topic  of  further  works. 
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